结合无人机载LiDAR点云法向量的聚类精简
结合无人机载LiDAR点云法向量的K-means++聚类精简
李沛婷1,2,3, 赵庆展
1.石河子大学信息科学与技术学院,石河子 832003
2.国家遥感中心新疆兵团分部,石河子 832003
3.兵团空间信息工程技术研究中心,石河子 832003
4.石河子大学机械电气工程学院,石河子 832003
摘要
点云精简可有效降低无人机载LiDAR数据量,对后期点云存储和快速处理具有重要意义。采用K-means++方法对点云法向量进行聚类,以实现点云精简。首先,利用回波次数去除多次回波点云,在使用零-均值标准化方法对点云属性归一化后,利用KD树(K-dimension tree)建立点云索引构建点云K邻域; 然后,采用主成分分析法估算点云法向量,借助肘方法确定最佳聚类数目; 最终,通过K-means++聚类方法实现点云精简。将精简结果生成Delaunay三角网并转换为栅格数据,通过相关系数验证方法的有效性。结果表明: 针对研究区69 544个点云数据,该方法可去除多次回波点云7 722个; 对点云法向量进行聚类数目为8的K-means++聚类,对应的精简率为分别为81.389%,81.833%和85.369%时效果较优; 精简后生成Delaunay三角网的时间远低于精简前,且当按81.833%进行精简处理时,相关系数最高,为0.890。该方法可为点云精简提供参考。
关键词： 点云K邻域 ; 点云法向量 ; K-means++聚类 ; Delaunay三角网
本文引用格式
李沛婷, 赵庆展, 田文忠, 马永建. 结合无人机载LiDAR点云法向量的K-means++聚类精简. 国土资源遥感[J], 2020, 32(2): 103-110 doi:10.6046/gtzyyg.2020.02.14
LI Peiting, ZHAO Qingzhan, TIAN Wenzhong, MA Yongjian. Point cloud simplification method combining K-means++ clustering with UAV LiDAR point cloud normal vectors. REMOTE SENSING FOR LAND & RESOURCES[J], 2020, 32(2): 103-110 doi:10.6046/gtzyyg.2020.02.14
0 引言无人机载激光雷达(light detection and ranging,LiDAR)是一种新型主动获取地表信息的观测系统,其集成激光扫描仪、惯性导航系统及计算机等技术[1],能够较好地克服传统方式获取地表信息的局限性,快速获取高精度地物数据,广泛应用于地形勘测、灾害监测、城市规划等领域[2]。但随着点云精度的提高,点云数据量也逐渐增大,这不仅增加数据存储负担,而且大大降低处理速度,故在保证不同应用精度的前提下,需对点云进行精简处理[3]。点云精简即以最少的数据量来表达地物的必要信息,是点云数据处理的基础[4]。
目前点云精简方法按照是否构建点云三角网分为直接采用点云和基于三角网格2种,基于三角网格的点云精简方法首先建立点云拓扑网格关系,然后合并细节较小的网格,去除对应点云,达到精简目的,主要包括包络网格法、顶点聚类法、区域合并法、边折叠法、小波分解等方法[5]。但存在建立点云拓扑结构速度慢的问题,且需按照一定阈值识别三角网以实现点云精简。直接采用点云的精简方法是根据空间位置关系计算出对应的离散几何信息实现精简[6],主要包括根据比例压缩、基于聚类压缩、基于曲率压缩、基于信息量压缩和基于网格等方法实现点云精简。葛源坤[7]针对平坦研究区采用包围盒法,针对突变区域采用基于点云曲率的最小距离法,结果表明将2种常用的点云精简方法结合可保证点云的关键几何信息; 陈西江等[8]首先利用主成分分析(principal component analysis,PCA)方法估计点云法向量及其夹角,然后利用最邻近点获取法向量夹角的信息熵,根据熵的大小实现点云精简; Han等[9]首先计算点云法向量,然后基于法向量识别和保留边缘点,删除非边缘点,实现点云精简; 苏本跃等[10]采用K-means对点云曲率进行聚类可以实现点云的去噪和精简处理; 唐林[11]提出保留点云属性的点云精简方法,即对所有的点云数据X-Y边界快速提取并保留,且对提取过边界的散乱点构建K邻域,通过邻域计算模糊熵和曲率,将曲率小于所设阈值的点剔除以实现点云精简。结合上述可知,实现点云精简方法是多样的,主要通过不同方法对点云的特征信息进行处理,以实现精简。其中K-means聚类方法可以实现点云精简处理,但仅对点云曲率进行聚类,对其他特征信息的聚类是否可实现聚类需进一步探讨。根据上述分析,本文在计算点云法向量的基础上,利用K-means++聚类方法对点云法向量聚类实现点云精简,并通过计算相关系数验证点云精简的有效性。
1 研究区概况及数据源1.1 研究区概况
以新疆第八师一五零团五连周边荒漠(E85°58'35″~85°59'4″,N44°57'29″~44°58'0″)作为无人机LiDAR数据获取的研究区,该荒漠区整体高程变化小,主要以旱生和沙生灌木为主,如梭梭、骆驼蓬、麻黄、碱蓬草和驼绒藜等。
以Scout B1-100单旋翼油动无人直升机作为飞行平台,搭载Riegl VUX-1激光扫描仪系统、OxTS Survey+2惯性导航系统和可见光相机,于2017年7月31日在研究区正上方以60 m的航高、6 m/s的巡航速度,550 kHz的扫描频率采集数据。激光扫描仪的基本参数见表1。Survey+2系统定位精度可达2 cm,横滚和俯仰精度为0.03°。
表1 激光扫描仪Riegl VUX-1主要参数
Tab.1 Parameters of Riegl laser scanner VUX-1
| 扫描仪参数 | 规格 |
|---|---|
| 扫描速度/(转·s-1) | 10~200 |
| 最小测量距离/m | 5 |
| 激光脉冲频率/kHz | 550 |
| 测量精度/mm | 15 |
| 角度分辨率/(°) | 0.001 |
| 最大视场角/(°) | 330 |
| 回波强度/bit | 16 |
| 回波次数 | 无限 |
在获取研究区扫描仪、惯性导航数据和可见光数据后,使用Riegl 公司提供的配套软件对点云数据进行航带裁剪拼接、扫描仪数据姿态矫正、点云去噪声点的预处理。截取部分实验区域作为本文的研究区,且以WGS_1984_UTM_Zone_44 N为投影坐标系导出点云数据,以.las格式存储点云。
1.2 数据源
经统计,研究区点云总个数为69 544,令X,Y,Z,In,Re,An分别表示点云三维坐标、回波强度值、回波次数和扫描角度,则研究区X范围为[-46.065 0,-15.002 3],Y范围为[83.958 8,112.144 5],Z范围为[275.863 0,280.533 3],In范围为[12 036,49 850],激光扫描仪系统记录的最大Re为5,An范围为[41°,54°]。研究区点云以回波强度进行三维可视化的结果如图1所示。
图1 研究区原始点云回波强度渲染结果
Fig.1 Visualization of the original point cloud echo intensity in study area
2 研究方法本研究技术流程如图2所示。
图2 本文技术流程
Fig.2 Technical flow chart for this paper
采用MATLAB R2017a软件提取X,Y,Z,In和Re,对X,Y,Z和In进行标准化以减少量纲引起的误差。在利用KD树(K-dimension tree)实现散乱点云的结构存储基础上,构建点云的K邻域,且利用PCA法估算得到点云法向量[12]。利用Python编程实现肘方法以确定K-means++方法的最佳聚类数目,进而采用K-means++方法对法向量进行聚类。接着,三维可视化不同聚类结果,通过目视解译获取保存研究区地物特征较好的簇,对簇中点云生成Delaunay三角网,对比参与聚类前点云生成的Delaunay三角网,且将三角网转换为栅格数据计算它们的相关系数,分析点云精简结果。
2.1 基于回波次数的点云预处理
点云回波分为单次回波和多次回波,其中单次回波多发生在地面或接近树冠的部分; 在多次回波中,首次回波大多为接近树冠部分的反射信号,末次回波一般是较树冠低一些的枝叶或地面的反射信号,除了首次和末次回波外的回波信号为中间次数的回波[13]。单次回波的回波次数为1,回波序号为1/1; 多次回波的回波序号是接收对应激光脉冲返回的顺序[14]。为了减少点云数据量,提高后期聚类速度,可以先利用点云的回波次数对点云进行预处理,去除多次回波点云。
聚类方法原理是根据不同距离来度量数据相似性,因此为了消除数据量纲和取值范围差异的影响,需对数据进行标准化[15],本文采用零-均值标准化方法将数据变为标准正态分布,标准化方程为
式中: μt和σt分别为点云属性t的均值和标准差; ti表示第i点的属性,包括点云三维坐标X,Y,Z和强度In; N为点云个数; χt和χ't分别表示属性标准前后的数据。根据上述分析,仅需对第1次回波点云中的三维坐标和回波强度值进行标准化处理。
2.2 点云K邻域构建
点云数据具有海量、离散、不规则等特点,为提高点云的搜索效率,需建立点云的数据结构。常用的数据结构包括八叉树、规则格网和KD树,其中KD树可以灵活地实现点云的空间索引,为构建邻域的提供索引[16]。KD树是在K维空间对数据的划分,把整个空间划分为特定的部分,并建立多维数据的索引,其每个结点对应K维空间的一个点[17],树的每层根据该层的分辨器做出分枝决策。其中,KD树第i层的分辨器定义为i除以K的余数,即imodK。KD树要么是一颗空树,要么是一颗满足条件的二叉树[18]。
2.3 点云法向量估算
对比最小二乘方法和PCA法估算法向量的结果,可知在平滑研究区PCA可以得到较高的精度[19],且结合本文研究区,本文选择PCA法估算点云法向量。其估算原理是从最小二乘法中推导得到,即根据点云中每个点p及其K邻域和最小二乘方法拟合点p的局部邻域得到局部平面β,平面β的数学表达式为[20]
式中: η为平面β的法向量; d为β到坐标原点的距离; K为点p最近K个点; argmin为当函数取最小值时对应自变量的取值;
2.4 K-means++聚类方法
2.4.1 K-means++聚类原理
通过K-means聚类可以较好地对点云进行精简[21],且K-means++聚类是对K-means聚类方法的改进,K-means++聚类通过使初始聚类中心点彼此尽可能远离的方法,使聚类结果更好、更一致,其初始聚类中心的选择步骤为:
1)初始化空数据集合M,存放k个聚类中心;
2)从输入数据中随机选择第一个中心点ci,将其加入集合M中;
3)对于集合M外的每一个样本点pi,计算pi到ci的距离D(pi),并找到D2(pi)最小的样本点p'i;
4)使用加权概率分布
5)重复步骤2)—4),至选定k个初始中心。
根据上述方法得到k个初始聚类中心ci,i=1,2,…,k。接着执行K-means算法,先设置最大迭代次数m和收敛的阈值ε。同时根据距离标准,计算每个样本点到中心ci之间的距离,且将该样本点划分到距离它最近的中心点ci,i∈1,2,…,k; 计算每个类别中相邻2次聚类中心变量D(pi),判断D(pi)是否小于收敛阈值ε,如果大于ε,且未达到最大迭代次数,则继续聚类; 否则终止聚类。
2.4.2 肘方法确定最佳聚类数目
K-means聚类方法可有效地对大数据进行伸缩和高效率处理,但是其难点之一是必须事先指定聚类数目,k值的选择会严重影响聚类结果[22],故本文使用肘方法确定最佳k值。基本思想是簇内误差平方和(within-cluster sum of squared errors,SSE)(即聚类偏差)骤增时的k值[23],即为最佳k值。即先对不同k值进行K-means聚类,根据聚类偏差得到不同k值的曲线,且绘制不同k值对应的聚类偏差图,最后找到聚类偏差骤增或者最显著拐点时对应的k值[24],其中聚类偏差SSE的定义为
式中: N为样本个数; k为聚类数目; pi为第i个样本点; cj为簇j的中心; 如果样本pi属于簇j,则w(i,j)=1,否则w(i,j)=0。
K-means聚类方法的另一个难点是初始聚类中心的选择。为了得到最优的聚类结果,本文在随机选择初始聚类中心和采用肘方法最佳聚类数目的基础上,针对不同的随机初始中心独立运行K-means聚类方法100次,接着从中选择SSE最小的作为最终模型。
2.5 点云精简验证
不规则三角网(triangulated irregular network,TIN)是点云的组织方式之一,其将点云按照一定区域的点云密度和位置划分为不规则的三角形网格。常见的TIN是Delaunay三角网,其满足最大最小角特性、最大空外接圆特性和唯一性[25]。最大最小角特性是指由点集形成的所有三角网中,Delaunay三角网中三角形的最小角度却是最大的。最大空外接圆特性是指由点集所形成的每一个三角形的外接圆均不包含点集中的其他任意点,即每个Delaunay三角网都是选择最近的点构建三角形。因为建立Delaunay三角网必须遵从以上2个特性,所以针对同一数据,不管采用什么建立三角网方法,其结果是唯一的。本文选择Delaunay三角网对点云进行重新组织可视化,且利用ArcGIS10.3软件将三角网转换为栅格数据,采用band collection statistics工具计算聚类前和聚类后对应栅格的相关系数。
3 结果与分析3.1 点云法向量估算结果
3.1.1 点云K邻域构建
研究区原始点云回波次数三维可视化如图3,通过分析图可知,第1次回波点云个数为61 822,第2—5次回波分别对应的点云个数为7 359,354,7和2,且第1次回波对应的点云包含了研究区全部地物,所以可去除多次回波点云,只针对第1次回波点云进行K-means++聚类实现点云精简。针对第1次回波对应的61 822个点云构建每个点云的K邻域,以欧氏距离作为度量标准,为了验证本文方法可行性,选取K=10以构建点云的K邻域,即每个点的邻域是由与该点距离最近的10个点组成。借助Python中Scikit-learn库,实现基于KD树来构建点云K邻域。为详细说明构建邻域的过程,本文选择部分点云进行说明,表2是索引号为1,2,3对应的K邻域,如以点云索引为1为例,即索引号为1的K邻域是由点云索引号为45,1,44,46,142,141,461,222,346和462构成的,其他点云类似。
图3 研究区原始点云回波次数三维可视化结果
Fig.3 3D visualization results of original point cloud echo times in study area
表2 部分点云K邻域对应的索引号
Tab.2 Index numbers corresponding to point cloud K neighborhoods
利用MATLAB软件编程找到每个点云邻域对应的三维坐标和回波强度值信息,表3是索引号为1对应K邻域的点云属性信息,其中Xstan,Ystan,Zstan和Instan分别为三维坐标和强度标准化后的数值。可以通过编程获取所有点云K邻域对应的属性信息,以实现点云K邻域的构建。
表3 索引号1对应的点云属性信息
Tab.3 Point cloud characteristic information corresponding to index number 1
| 索引号 | Xstan | Ystan | Zstan | Instan |
|---|---|---|---|---|
| 45 | -1.362 | -1.096 | 0.308 | 1.097 |
| 1 | -1.340 | -1.101 | 0.304 | 0.968 |
| 44 | -1.378 | -1.085 | 0.329 | 0.983 |
| 46 | -1.344 | -1.108 | 0.323 | 1.198 |
| 142 | -1.345 | -1.113 | 0.337 | 1.058 |
| 141 | -1.360 | -1.103 | 0.347 | 1.140 |
| 461 | -1.380 | -1.118 | 0.315 | 1.239 |
| 222 | -1.375 | -1.099 | 0.355 | 1.221 |
| 346 | -0.609 | -1.632 | 0.620 | 0.179 |
| 462 | -1.362 | -1.131 | 0.314 | 1.094 |
3.1.2 点云法向量估算结果
对点云邻域进行PCA分析,找到最小特征值对应的特征向量以实现点云法向量的估算,利用MATLAB软件编程实现点云法向量估算。为可视化点云法向量,本文仅显示部分区域的点云法向量,选取区域的X,Y和Z范围分别为[-41,-21],[87,106]和[276,280],法向量可视化结果如图4所示。
图4 点云法向量可视化
Fig.4 Visualization of point cloud normal vector
3.2 点云精简结果分析
3.2.1 K-means++聚类结果
基于不同随机初始中心独立运行K-means方法100次,每次运行最大迭代次数为3 000次,选择聚类偏差最小的作为最终聚类结果; 容忍度设为0.000 1,即判断是否收敛的阈值为0.000 1。首先采用Python语言编程实现肘方法以确定K-means++对点云法向量聚类的最佳聚类数目,图5为聚类偏差结果,横坐标表示聚类数目,纵坐标表示聚类偏差。分析图5可知,当k=8时,聚类偏差呈现肘型,即对于点云法向量属性,最佳聚类数目为8。
图5 聚类偏差结果
Fig.5 Clustering SSE for different number of results
K-means++对法向量进行聚类数目为8的聚类,独立运行K-means++方法100次需消耗的时间为22.392 s。K-means++对法向量聚类的三维可视化结果如图6所示,图中不同颜色表示点云不同强度值,图6(a)—(h)分别对应聚类结果中的簇1—簇8。
图6 K-means++对法向量聚类的三维可视化
Fig.6 3D visualization of K-means++ for clustering point cloud normal vector
对比聚类前后点云结果可知,簇1、簇2和簇3保持研究区地物特征,其他簇虽然很好地保存地形,但地物边界保留不完整。对聚类结果中点云个数和强度误差进行统计得到表4,分析表可知,因为簇2和簇3基本包含了研究区地物特征,所以对应的最小强度误差为0,且簇3强度误差范围最小,为87。结合聚类结果可视化和统计可知,按照精简率为81.389%,85.369%和81.833%进行点云精简处理比按照精简率为90.916%,88.300%,89.606%,94.031%,88.557%处理的效果更好,即点云越多更容易表达地物细节特征。
表4 聚类后点云个数和强度误差统计
Tab.4 Statistics of point cloud number and intensity error after clustering
| 类别 | 聚类后的 点云个数 | 点云精简 率/% | 点云回波强度误差 | |
|---|---|---|---|---|
| 最小误差 | 最大误差 | |||
| 簇1 | 11 506 | 81.389 | 437 | 787 |
| 簇2 | 9 045 | 85.369 | 0 | 961 |
| 簇3 | 11 231 | 81.833 | 0 | 87 |
| 簇4 | 5 616 | 90.916 | 677 | 2 534 |
| 簇5 | 7 234 | 88.300 | 437 | 612 |
| 簇6 | 6 426 | 89.606 | 524 | 2 469 |
| 簇7 | 3 690 | 94.031 | 1 224 | 1 551 |
| 簇8 | 7 074 | 88.557 | 852 | 918 |
3.2.2 点云精简分析
通过上述分析,将簇1、簇2、簇3和聚类前对应的点云生成Delaunay三角网。接着利用自然邻域插值将三角网转换为栅格数据,其中栅格像元大小均设置为0.125,且分别计算簇1、簇2和簇3与聚类前点云栅格数据的相关系数 。图7为Delaunay三角网可视化结果。统计聚类前后Delaunay三角网,并计算对应的相关系数,结果见表5。
图7 簇1、簇2和簇3对应点云的Delaunay三角网
Fig.7 Delaunay triangulation of clusters 1, clusters 2 and clusters 3
表5 聚类前后Delaunay三角网统计结果和相关系数
Tab.5 Statistical results and correlation coefficients of Delaunay triangulation before and after clustering
| 点云精简 率/% | 类别 | 三角网 数据/个 | 三角网顶 点数据/个 | 构建三角 网时间/s | 相关系数 |
|---|---|---|---|---|---|
| 81.389 | 簇1 | 22 983 | 11 505 | 0.28 | 0.884 |
| 85.369 | 簇2 | 18 061 | 9 044 | 0.32 | 0.870 |
| 81.833 | 簇3 | 22 434 | 11 230 | 0.23 | 0.890 |
| 100 | 聚类前 | 123 578 | 61 809 | 1.01 | 1 |
分析表5可知,采用精简后点云生成三角网时间大大缩短; 按照81.389%,85.369%和81.833%进行精简处理后的点云基本参与三角网的生成中; 当点云精简率为81.833%时,其对应的相关系数最高0.890,且运行时间最短,为0.23 s。
1)利用肘方法可以较好地确定聚类数目; 同时,随机抽取点云作为初始中心,且独立运行聚类方法多次,结合聚类偏差定量选取最终聚类结果,较好地解决聚类两大问题。
2)通过KD树建立点云的数据结构,利用其构建点云的K邻域,且采用PCA方法估算法向量。最后利用K-means++聚类方法实现点云精简,精简率分别为81.389%,85.369%和81.833%。
3)通过生成Delaunay三角网,且转换为栅格数据,计算精简前后的相关系数。综合考虑运行时间和相关系数,可知点云精简率为81.833%时效果最好。
4)本文方法通过一次点云聚类处理即可获取不同精简率的点云,根据点云应用领域的不同可以随机选取不同的精简结果。本文方法的不足之处在于: 只关注K-means++聚类方法可以实现点云的精简处理,且验证精简结果,但在建立点云K邻域时,本文直接选择K=10,虽然可以提高邻域构建速度,但是对不同邻域K值的设定影响点云法向量的估算结果,需要进一步研究。
参考文献（略）